Aplicaciones de Agro Data Science

Análisis de mortalidad de insectos y predicción de rendimiento

Christian R.A. Vásquez Velasco

InkaStats Academy

2024-09-15

Agro Data Science

La agricultura moderna


Agro Data Science


  • ADS combina el conocimiento agronómico con las tecnologías de la información para resolver problemas en la agricultura.
  • Abarca temas desde la recopilación y análisis de datos hasta la visualización y comunicación de resultados.
  • Los científicos de datos agrícolas trabajan con agricultores, empresas agroindustriales, universidades y organizaciones gubernamentales para mejorar la eficiencia y la sostenibilidad de la agricultura.

Problemas que resuelven el ADS


  • Gestión del riesgo climático
  • Determinación de dosis y momento de fertilización
  • Análisis y predicción de rendimiento
  • Estrés hídrico
  • Monitoreo de campos
  • Diagnóstico de plagas y enfermedades

  • Importancia del acceso a datos precisos y confiables: Estaciones meteorológicas, sensores de campo, drones y satélites.
  • Herramientas y técnicas de análisis de datos: Visualización de datos, Clasificación y Regresión, Series temporales, Aprendizaje automático, Inteligencia artificial.

Desafios

  • Falta de acceso a datos de calidad.
  • Cuestiones éticas y de privacidad.

Oportunidades

  • Mejorar la eficiencia y la sostenibilidad de la agricultura.
  • Tomar decisiones informadas y basadas en datos.

R

  • R es un lenguaje de programación utilizado en el análisis estadístico y la visualización de datos.
  • Fue creado en 1993 por Ross Ihaka y Robert Gentleman en la Universidad de Auckland, Nueva Zelanda.
  • R es un software libre y de código abierto.
  • R cuenta con una gran cantidad de paquetes que se utilizan para realizar tareas específicas en el análisis de datos.
  • R tiene una sintaxis fácil de aprender y es utilizado tanto por científicos de datos como por estadísticos, ingenieros y profesionales de negocios.

Tidyverse

  • Tidyverse es un conjunto de paquetes de software de análisis de datos en el lenguaje de programación R, que se utilizan para manipular, visualizar y modelar datos.
  • Estos paquetes comparten una filosofía común y una sintaxis consistente, lo que hace que sea fácil trabajar con ellos en conjunto.
  • Los paquetes que forman parte de Tidyverse incluyen dplyr, ggplot2, tidyr, readr, purrr, tibble, stringr, forcats, entre otros.

Tidymodels

  • Tidymodels es un conjunto de paquetes de software de modelado estadístico en el lenguaje de programación R que se utilizan para construir, ajustar, validar y analizar modelos predictivos.
  • El objetivo principal de Tidymodels es proporcionar herramientas para ayudar a los usuarios a realizar modelado estadístico de manera más eficiente y efectiva, utilizando técnicas modernas de modelado y validación..
  • Los paquetes que forman parte de Tidymodels incluyen recipes, parsnip, yardstick, dials, infer, entre otros.

Análisis de mortalidad de insectos

Entendiendo el problema

Las áreas de sanidad cuentan con data histórica.

Deben transformarla en información útil.

¿?

Predecir la mortalidad que se obtendrá en las aplicaciones futuras.

Para ilustración se utilizará el inventario de Eficacia de plaguicidas.

Data summary
Name data
Number of rows 929
Number of columns 6
_______________________
Column type frequency:
factor 3
numeric 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
InsectStage 0 1 FALSE 134 133: 41, 134: 41, 62: 39, 61: 32
Method 0 1 FALSE 7 5: 707, 1: 133, 4: 75, 6: 6
Product 0 1 FALSE 101 1: 145, 94: 38, 64: 26, 73: 26

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Pre_Mean_Count 0 1.00 83.41 177.35 0 6.41 24.23 85.67 2416 ▇▁▁▁▁
Post_Mean_Count 5 0.99 32.05 123.41 0 0.22 2.33 18.19 1929 ▇▁▁▁▁
Mortality_Mean_Cor 0 1.00 0.63 0.41 0 0.15 0.84 0.98 1 ▃▁▁▂▇

Método: método de aplicación
- 1 Control - 2 Insecticidas en gránulos - 3 Tratamiento de semillas - 4 Aplicación al suelo - 5 Pulverización - 6 Niebla de vapor - 7 Aplicación superficial con agua

# A tibble: 10 × 6
   InsectStage Method Product Pre_Mean_Count Post_Mean_Count Mortality_Mean_Cor
   <fct>       <fct>  <fct>            <dbl>           <dbl>              <dbl>
 1 1           1      1                87.8           150.                0    
 2 1           3      57                1.67            3                 0    
 3 1           3      71                2.67            2.67              0    
 4 1           4      31               45               6.33              0.859
 5 1           4      49               45               5.67              0.874
 6 1           4      57               90.6            48.6               0.464
 7 1           4      78               39.3            16.3               0.585
 8 1           5      58               74.3            46.7               0.372
 9 2           1      1                13.3            53.7               0    
10 2           4      57                8              11                 0    

Partición inicial

split <- initial_split(data, prop = 0.70, strata = Mortality_Mean_Cor)

Partición inicial

split <- initial_split(data, prop = 0.70, strata = Mortality_Mean_Cor)
data_train <- split %>% training()
tibble [649 × 6] (S3: tbl_df/tbl/data.frame)
 $ InsectStage       : Factor w/ 140 levels "1","2","3","4",..: 1 2 2 2 3 3 4 4 5 7 ...
 $ Method            : Factor w/ 7 levels "1","2","3","4",..: 1 1 4 5 1 3 1 4 1 1 ...
 $ Product           : Factor w/ 101 levels "1","2","3","4",..: 1 1 57 58 1 57 1 58 1 1 ...
 $ Pre_Mean_Count    : num [1:649] 87.8 13.3 8 14.7 117.3 ...
 $ Post_Mean_Count   : num [1:649] 150.5 53.7 11 18.3 167.9 ...
 $ Mortality_Mean_Cor: num [1:649] 0 0 0 0 0 ...

Partición inicial

split <- initial_split(data, prop = 0.70, strata = Mortality_Mean_Cor)
data_train <- split %>% training()
data_test <- split %>% testing()
tibble [280 × 6] (S3: tbl_df/tbl/data.frame)
 $ InsectStage       : Factor w/ 140 levels "1","2","3","4",..: 1 1 1 1 1 3 4 6 6 6 ...
 $ Method            : Factor w/ 7 levels "1","2","3","4",..: 3 3 4 4 4 3 4 1 5 5 ...
 $ Product           : Factor w/ 101 levels "1","2","3","4",..: 57 71 49 57 78 71 57 1 14 69 ...
 $ Pre_Mean_Count    : num [1:280] 1.67 2.67 45 90.6 39.33 ...
 $ Post_Mean_Count   : num [1:280] 3 2.67 5.67 48.6 16.33 ...
 $ Mortality_Mean_Cor: num [1:280] 0 0 0.874 0.464 0.585 ...

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train)

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`))

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("Mortality_Mean_Cor")),
                 fn = ~log(.))

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("Mortality_Mean_Cor")),
                 fn = ~log(.))

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("Mortality_Mean_Cor")),
                 fn = ~log(.)) %>%
  step_dummy(Method, one_hot = TRUE) %>%
  step_dummy(Product, one_hot = TRUE) %>%
  step_dummy(InsectStage, one_hot = TRUE)

Receta de preprocesamiento

data_prep <- prep(data_rec)

Receta de preprocesamiento

data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)

Receta de preprocesamiento

data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)
data_test_prep  <- bake(data_prep, new_data = data_test)

Receta de preprocesamiento

data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)
data_test_prep  <- bake(data_prep, new_data = data_test)
head(data_train_prep)
# A tibble: 6 × 250
  Pre_Mean_Count Mortality_Mean_Cor Method_X1 Method_X2 Method_X3 Method_X4
           <dbl>              <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
1           4.48                  0         1         0         0         0
2           2.59                  0         1         0         0         0
3           2.08                  0         0         0         0         1
4           2.69                  0         0         0         0         0
5           4.76                  0         1         0         0         0
6           2.46                  0         0         0         1         0
# ℹ 244 more variables: Method_X5 <dbl>, Method_X6 <dbl>, Method_X7 <dbl>,
#   Product_X1 <dbl>, Product_X2 <dbl>, Product_X3 <dbl>, Product_X4 <dbl>,
#   Product_X5 <dbl>, Product_X6 <dbl>, Product_X7 <dbl>, Product_X8 <dbl>,
#   Product_X9 <dbl>, Product_X10 <dbl>, Product_X11 <dbl>, Product_X12 <dbl>,
#   Product_X13 <dbl>, Product_X14 <dbl>, Product_X15 <dbl>, Product_X16 <dbl>,
#   Product_X17 <dbl>, Product_X18 <dbl>, Product_X19 <dbl>, Product_X20 <dbl>,
#   Product_X21 <dbl>, Product_X22 <dbl>, Product_X23 <dbl>, …

Random Forest

Definir especificaciones del modelo

# Inicializar un objeto de random forest
rf_model_spec <- rand_forest(
  trees = tune(), 
  min_n = tune(), 
  mtry = tune()
) %>% 
  # Establecer el modelo de motor
  set_engine('ranger',                         
             importance = "permutation") %>% 
  # Establecer el modo de modelo
  set_mode('regression')

Definir fórmula del modelo

formula_rf = Mortality_Mean_Cor ~ .  

Definir flujo de trabajo

rf_wf <- workflow() %>%
  add_recipe(data_rec) %>%
  add_model(rf_model_spec,
            formula = formula_rf) 
#obtener todos los argumentos ajustables posibles en el flujo de trabajo
tune_args(rf_wf) 
# A tibble: 3 × 6
  name  tunable id    source     component   component_id
  <chr> <lgl>   <chr> <chr>      <chr>       <chr>       
1 mtry  TRUE    mtry  model_spec rand_forest <NA>        
2 trees TRUE    trees model_spec rand_forest <NA>        
3 min_n TRUE    min_n model_spec rand_forest <NA>        

Extraer el conjunto de hiperparámetros del modelo

rf_params <- hardhat::extract_parameter_set_dials(rf_wf) %>% 
  update(trees = trees(range = c(50L, 150L)))
rf_params <- rf_params %>% 
  dials::finalize(rf_params)
rf_params
Collection of 3 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[+]
      trees trees nparam[+]
      min_n min_n nparam[+]

Definir métricas del performance predictivo

multi_met <- yardstick::metric_set(
yardstick::rmse,
yardstick::rsq,
yardstick::mape,
yardstick::mae)

Definir esquema de validación cruzada

set.seed(234)
data_folds <- rsample::vfold_cv(data_train,
                                 v = 20)

Búsqueda del mejor conjunto de Hiperparámtros mediante validación cruzada

set.seed(2020)
tictoc::tic()
rf_tune <-
  rf_wf %>%
  tune_bayes(
    resamples = data_folds,
    param_info = rf_params,
    # por defecto es 5 (número aleatorio de combinaciones (puntos de grilla) de hiperparámetros)
    initial = 5, 
    iter = 50,
    metrics = multi_met,
    control = control_bayes(
      # El corte entero para el número de iteraciones sin mejores resultados.
      no_improve = 10,
      extract = identity,
      save_pred = TRUE,
      verbose = FALSE
      )
  )
tictoc::toc()
320 sec elapsed

Métricas del rendimiento predictivo en fase de validación cruzada

Seleccionar el mejor conjunto de hiperparámetros

best_rf_model <- select_best(rf_tune, "rmse")
final_rf <- finalize_workflow(rf_wf,
                              best_rf_model)
final_rf
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_select()
• step_mutate_at()
• step_dummy()
• step_dummy()
• step_dummy()

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (regression)

Main Arguments:
  mtry = 6
  trees = 76
  min_n = 39

Engine-Specific Arguments:
  importance = permutation

Computational engine: ranger 

Métricas del rendimiento predictivo para el mejor conjunto de hiperparámetros

# Métricas promedio de todas las particiones
metrics_rf <- 
  rf_tune %>% 
  collect_metrics(summarize = TRUE) %>%
  dplyr::mutate(modelo = "rf") %>%
  dplyr::select(-c(.estimator,.config)) %>%
  dplyr::filter(mtry %in% best_rf_model$mtry,
                min_n %in% best_rf_model$min_n,
                trees %in% best_rf_model$trees)
metrics_rf
mtry trees min_n .metric mean n std_err .iter modelo
6 76 39 mae 0.2134 20 0.004847 11 rf
6 76 39 mape Inf 20 NaN 11 rf
6 76 39 rmse 0.2688 20 0.006324 11 rf
6 76 39 rsq 0.5735 20 0.024704 11 rf

Entrenamiento del modelo final

rf_model <- final_rf %>%
  fit(data_train) %>%
  extract_fit_parsnip()
rf_model
parsnip model object

Ranger result

Call:
 ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~6L,      x), num.trees = ~76L, min.node.size = min_rows(~39L, x),      importance = ~"permutation", num.threads = 1, verbose = FALSE,      seed = sample.int(10^5, 1)) 

Type:                             Regression 
Number of trees:                  76 
Sample size:                      649 
Number of independent variables:  249 
Mtry:                             6 
Target node size:                 39 
Variable importance mode:         permutation 
Splitrule:                        variance 
OOB prediction error (MSE):       0.07622 
R squared (OOB):                  0.5352 

Importancia de variables

Testeo interno

Resumen de la predicción

predicciones <- rf_model %>%
  predict(
    new_data = data_test_prep,
    type = "numeric"
  )
  
predicciones_rf <- predicciones %>% 
  bind_cols(data_test_prep %>% dplyr::select(Mortality_Mean_Cor)) %>%
  dplyr::mutate(modelo = "rf")

summary(predicciones_rf)  
     .pred       Mortality_Mean_Cor    modelo         
 Min.   :0.127   Min.   :0.000      Length:280        
 1st Qu.:0.530   1st Qu.:0.135      Class :character  
 Median :0.719   Median :0.847      Mode  :character  
 Mean   :0.635   Mean   :0.623                        
 3rd Qu.:0.809   3rd Qu.:0.983                        
 Max.   :0.919   Max.   :1.000                        

Testeo interno

Resumen de las métricas de predicción

error_test_rf  <- multi_met(
  data     = predicciones_rf,
  truth    = Mortality_Mean_Cor,
  estimate = .pred,
  na_rm    = TRUE
) %>%
  dplyr::mutate(
    modelo = "rf"
  )
error_test_rf
# A tibble: 4 × 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard       0.275 rf    
2 rsq     standard       0.591 rf    
3 mape    standard     Inf     rf    
4 mae     standard       0.225 rf    

Testeo interno

Predicción de rendimiento

Fórmulas para la proyección del rendimiento

Debemos entender que de forma empírica, las estimaciones de cosecha se realizan en muchos cultivos teniendo en cuenta una fórmula matemática como la siguiente:

\[Producción = A * DS * NOP * PO\]

Dónde:

A = Área

DS = Densidad de siembra

NOP = Número de órganos por planta

PO = Peso de órganos

Entendiéndose con esta fórmula que para estimar la producción por campaña de cualquier cultivo, se necesita información de la extensión de siembra, número de plantas por hectárea, número de órganos a cosechar y el peso promedio de estos órganos.

Para fines logísticos, hacer una estimación de producción a nivel de campaña no ayuda mucho, ya que se necesita saber la producción en distintos puntos de corte o periodos de tiempo.

Para esto último fin, podemos hacer una modificación a la ecuación, obteniendo lo siguiente:

\[Producción_{(t)} = A * DS * NOP_{(t)} * PO_{(t)} * TM_{(t)}\]

Dónde:

A = Área

DS = Densidad de siembra

NOP = Número de órganos por planta

PO = Peso de órganos

TM = Tasa de maduración

Esta última ecuación añade un modelamiento temporal, que está en función de la tasa de maduración y crecimiento de los órganos a cosechar.

Entonces, para realizar Crop Yield Forecast, debemos entender que existe en el tiempo uno o más factores que condicionarán la producción en distintos puntos transversales y por ello debemos tomar datos horizontales o de series de tiempo en el cultivo, pues serán fundamentales para realizar las proyecciones. Es así, que debemos imaginar cualquier algoritmo como una función similar a esta última ecuación matemática y debemos lograr que el modelo por crear, a través del aprendizaje estadístico, cumpla con la lógica necesaria para simular este comportamiento matemático (que fue pensado para simular el comportamiento real del cultivo).

Wild Blueberry Yield Prediction

El artículo “Wild blueberry yield prediction using a combination of computer simulation and machine learning algorithms”, usa datos recolectados por 30 años en Maine (EEUU) y publicados por Wild Blueberry Pollination Model, para un modelo de simulación espacial. El modelamiento predictivo en arándanos requiere de datos con fuerte influyencia de factores espaciales (geográficos), y adicionalmente en este caso, de las plantas, especies de abejas y meteorología. Estos datos simulados son usados por investigadores para probar la calidad predictiva de sus algoritmos de aprendizaje estadístico a datos reales y datos generados por simulación. Emplearemos este conjunto de datos para entrenar modelos en el enfoque tidy.

Descriptive Statistics  

                    andrena   AverageOfLowerTRange   AverageOfUpperTRange   AverageRainingDays
----------------- --------- ---------------------- ---------------------- --------------------
             Mean      0.47                  48.61                  68.72                 0.32
          Std.Dev      0.16                   5.42                   7.68                 0.17
              Min      0.00                  41.20                  58.20                 0.06
               Q1      0.38                  45.80                  64.70                 0.10
           Median      0.50                  50.80                  71.90                 0.26
               Q3      0.63                  50.80                  71.90                 0.39
              Max      0.75                  55.90                  79.00                 0.56
              MAD      0.18                   7.41                  10.53                 0.24
              IQR      0.25                   5.00                   7.20                 0.29
               CV      0.34                   0.11                   0.11                 0.54
         Skewness      0.19                  -0.01                  -0.02                 0.07
      SE.Skewness      0.09                   0.09                   0.09                 0.09
         Kurtosis     -0.67                  -1.34                  -1.35                -1.28
          N.Valid    777.00                 777.00                 777.00               777.00
        Pct.Valid    100.00                 100.00                 100.00               100.00

Table: Table continues below

 

                    bumbles   clonesize   fruitmass   fruitset   honeybee   MaxOfLowerTRange
----------------- --------- ----------- ----------- ---------- ---------- ------------------
             Mean      0.28       18.77        0.45       0.50       0.42              59.31
          Std.Dev      0.07        7.00        0.04       0.08       0.98               6.65
              Min      0.00       10.00        0.31       0.19       0.00              50.20
               Q1      0.25       12.50        0.42       0.45       0.25              55.80
           Median      0.25       12.50        0.45       0.51       0.25              62.00
               Q3      0.38       25.00        0.48       0.56       0.50              66.00
              Max      0.58       40.00        0.54       0.65      18.43              68.20
              MAD      0.00        0.00        0.04       0.08       0.00               9.19
              IQR      0.13       12.50        0.06       0.11       0.25              10.20
               CV      0.23        0.37        0.09       0.16       2.35               0.11
         Skewness      0.15        0.57       -0.10      -0.52      16.70              -0.02
      SE.Skewness      0.09        0.09        0.09       0.09       0.09               0.09
         Kurtosis      1.68       -0.64       -0.47       0.05     296.31              -1.35
          N.Valid    777.00      777.00      777.00     777.00     777.00             777.00
        Pct.Valid    100.00      100.00      100.00     100.00     100.00             100.00

Table: Table continues below

 

                    MaxOfUpperTRange   MinOfLowerTRange   MinOfUpperTRange    osmia   RainingDays
----------------- ------------------ ------------------ ------------------ -------- -------------
             Mean              82.28              28.69              49.70     0.56         18.31
          Std.Dev               9.19               3.21               5.60     0.17         12.12
              Min              69.70              24.30              39.00     0.00          1.00
               Q1              77.40              27.00              46.80     0.50          3.77
           Median              86.00              30.00              52.00     0.63         16.00
               Q3              89.00              30.00              52.00     0.75         24.00
              Max              94.60              33.00              57.20     0.75         34.00
              MAD              12.75               4.45               7.71     0.18         18.13
              IQR              11.60               3.00               5.20     0.25         20.23
               CV               0.11               0.11               0.11     0.30          0.66
         Skewness              -0.01              -0.01              -0.02    -0.92         -0.21
      SE.Skewness               0.09               0.09               0.09     0.09          0.09
         Kurtosis              -1.35              -1.35              -1.34     0.54         -1.24
          N.Valid             777.00             777.00             777.00   777.00        777.00
        Pct.Valid             100.00             100.00             100.00   100.00        100.00

Table: Table continues below

 

                      Row#    seeds     yield
----------------- -------- -------- ---------
             Mean   388.00    36.12   6012.85
          Std.Dev   224.44     4.38   1356.96
              Min     0.00    22.08   1637.70
               Q1   194.00    33.12   5124.85
           Median   388.00    36.17   6107.38
               Q3   582.00    39.24   7022.19
              Max   776.00    46.59   8969.40
              MAD   287.62     4.56   1417.36
              IQR   388.00     6.12   1897.33
               CV     0.58     0.12      0.23
         Skewness     0.00    -0.06     -0.32
      SE.Skewness     0.09     0.09      0.09
         Kurtosis    -1.20    -0.40     -0.39
          N.Valid   777.00   777.00    777.00
        Pct.Valid   100.00   100.00    100.00

Partición inicial

split <- initial_split(data, prop = 0.70, strata = yield)

Partición inicial

split <- initial_split(data, prop = 0.70, strata = yield)
data_train <- split %>% training()
'data.frame':   541 obs. of  18 variables:
 $ Row#                : int  1 3 4 5 6 12 13 15 16 17 ...
 $ clonesize           : num  37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 ...
 $ honeybee            : num  0.75 0.75 0.75 0.75 0.75 0.25 0.25 0.25 0.25 0.25 ...
 $ bumbles             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ andrena             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ osmia               : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ MaxOfUpperTRange    : num  86 94.6 86 86 94.6 86 94.6 86 86 94.6 ...
 $ MinOfUpperTRange    : num  52 57.2 52 52 57.2 52 57.2 52 52 57.2 ...
 $ AverageOfUpperTRange: num  71.9 79 71.9 71.9 79 71.9 79 71.9 71.9 79 ...
 $ MaxOfLowerTRange    : num  62 68.2 62 62 68.2 62 68.2 62 62 68.2 ...
 $ MinOfLowerTRange    : num  30 33 30 30 33 30 33 30 30 33 ...
 $ AverageOfLowerTRange: num  50.8 55.9 50.8 50.8 55.9 50.8 55.9 50.8 50.8 55.9 ...
 $ RainingDays         : num  1 1 24 34 24 1 16 24 34 24 ...
 $ AverageRainingDays  : num  0.1 0.1 0.39 0.56 0.39 0.1 0.26 0.39 0.56 0.39 ...
 $ fruitset            : num  0.444 0.408 0.354 0.31 0.284 ...
 $ fruitmass           : num  0.425 0.409 0.383 0.366 0.352 ...
 $ seeds               : num  33.4 31.6 28.9 27.3 26.1 ...
 $ yield               : num  4948 4304 3436 2825 2625 ...

Partición inicial

split <- initial_split(data, prop = 0.70, strata = yield)
data_train <- split %>% training()
data_test <- split %>% testing()
'data.frame':   236 obs. of  18 variables:
 $ Row#                : int  0 2 7 8 9 10 11 14 19 27 ...
 $ clonesize           : num  37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 25 ...
 $ honeybee            : num  0.75 0.75 0.75 0.75 0.75 0.75 0.25 0.25 0.25 0.25 ...
 $ bumbles             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ andrena             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ osmia               : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ MaxOfUpperTRange    : num  86 94.6 94.6 77.4 77.4 69.7 86 94.6 77.4 94.6 ...
 $ MinOfUpperTRange    : num  52 57.2 57.2 46.8 46.8 42.1 52 57.2 46.8 57.2 ...
 $ AverageOfUpperTRange: num  71.9 79 79 64.7 64.7 58.2 71.9 79 64.7 79 ...
 $ MaxOfLowerTRange    : num  62 68.2 68.2 55.8 55.8 50.2 62 68.2 55.8 68.2 ...
 $ MinOfLowerTRange    : num  30 33 33 27 27 24.3 30 33 27 33 ...
 $ AverageOfLowerTRange: num  50.8 55.9 55.9 45.8 45.8 41.2 50.8 55.9 45.8 55.9 ...
 $ RainingDays         : num  16 16 34 16 1 16 16 1 16 16 ...
 $ AverageRainingDays  : num  0.26 0.26 0.56 0.26 0.1 0.26 0.26 0.1 0.26 0.26 ...
 $ fruitset            : num  0.411 0.384 0.247 0.428 0.464 ...
 $ fruitmass           : num  0.408 0.399 0.343 0.415 0.436 ...
 $ seeds               : num  31.7 30.5 25 32.3 34.8 ...
 $ yield               : num  3813 3867 2380 4235 5357 ...

Receta de preprocesamiento

data_rec <- recipe(formula = `yield` ~ .,
                   x = data_train) %>%
  step_select(where(base::is.numeric)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("RainingDays",
                             "AverageRainingDays",
                             "honeybee",
                             "bumbles",
                             "andrena",
                             "osmia",
                             "fruitset",
                             "fruitmass")),
                 fn = ~log(.))
                 
data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)
data_test_prep  <- bake(data_prep, new_data = data_test)

Definir función para recuperar coeficientes

get_glm_coefs <- function(x) {
  x %>% 
    # get the glm model object
    extract_fit_engine() %>% 
    # transform its format
    broom::tidy()
}

Definir validación cruzada

set.seed(234)
data_folds <- rsample::vfold_cv(data_train,
                                 v = 20)

Regresión múltiple

Definir especificaciones del modelo

# Inicializar un objeto de regresión lineal
linear_model_spec <- linear_reg() %>% 
  # Establecer el modelo de motor
  set_engine('lm') %>% 
  # Establecer el modo de modelo
  set_mode('regression')

Definir fórmula del modelo

formula_lm = yield ~ fruitmass + fruitset + seeds + osmia + andrena + bumbles +
  MinOfUpperTRange + clonesize + AverageRainingDays + honeybee + RainingDays 

Definir flujo de trabajo

lm_wf <- workflow() %>%
  add_recipe(data_rec) %>%
  add_model(linear_model_spec,
            formula = formula_lm) 

Validación cruzada

tictoc::tic()
validacion_lm <-
  lm_wf %>%
  fit_resamples(
    resamples = data_folds,
    metrics = multi_met,
    control = control_resamples( 
      extract = get_glm_coefs,
      save_pred = TRUE,
      verbose = FALSE
    )
  )
tictoc::toc()
13.86 sec elapsed

Métricas del rendimiento predictivo en fase de validación cruzada

# Métricas promedio de todas las particiones
metrics_lm <- 
  validacion_lm %>%
  unnest(.predictions) %>%
  group_by(id) %>%
  multi_met(truth = `yield`,
            estimate = `.pred`) %>%
  group_by(.metric) %>%
  dplyr::summarise(mean = mean(`.estimate`, na.rm = T),
            std_err = sd(`.estimate`, na.rm = T),
            n = n()) %>%
  ungroup() %>%
  dplyr::mutate(modelo = "lm")
metrics_lm
.metric mean std_err n modelo
mae 119.4882 15.619637 20 lm
mape 2.1658 0.385811 20 lm
rmse 151.1992 27.588467 20 lm
rsq 0.9881 0.005662 20 lm

Entrenamiento del modelo final

# Métricas promedio de todas las particiones
lm_model <- linear_model_spec %>%
  fit(formula_lm,
      data = data_train_prep)
lm_model
parsnip model object


Call:
stats::lm(formula = yield ~ fruitmass + fruitset + seeds + osmia + 
    andrena + bumbles + MinOfUpperTRange + clonesize + AverageRainingDays + 
    honeybee + RainingDays, data = data)

Coefficients:
       (Intercept)           fruitmass            fruitset               seeds  
          1.984826           -5.944940            2.612865            2.204394  
             osmia             andrena             bumbles    MinOfUpperTRange  
          0.070036            0.063956           -0.058212            0.002545  
         clonesize  AverageRainingDays            honeybee         RainingDays  
          0.031348           -0.035562            0.003223           -0.000328  
term estimate std.error statistic p.value
(Intercept) 1.9848262 0.4838516 4.1021 0.000047388210623994374997604295263187168529839255
fruitmass -5.9449399 0.3819019 -15.5667 0.000000000000000000000000000000000000000000002951
fruitset 2.6128649 0.1722453 15.1694 0.000000000000000000000000000000000000000000204388
seeds 2.2043935 0.1860721 11.8470 0.000000000000000000000000000069942530762398219600
osmia 0.0700361 0.0100628 6.9599 0.000000000010115665607954346622921262821570564938
andrena 0.0639560 0.0090381 7.0762 0.000000000004724348408114360337636122855187181813
bumbles -0.0582122 0.0264683 -2.1993 0.028287426461342017602884624238868127577006816864
MinOfUpperTRange 0.0025452 0.0194022 0.1312 0.895680124173958036593035103578586131334304809570
clonesize 0.0313481 0.0065414 4.7923 0.000002145610420601447544935308497926484960771631
AverageRainingDays -0.0355625 0.0603067 -0.5897 0.555647841162889788435563787061255425214767456055
honeybee 0.0032228 0.0026613 1.2110 0.226440579572288552823167151473171543329954147339
RainingDays -0.0003279 0.0008122 -0.4037 0.686618002087163170621408880833769217133522033691

Testeo interno

Resumen de la predicción

predicciones <- lm_model %>%
  predict(
    new_data = data_test_prep,
    type = "numeric"
  )
  
predicciones_lm <- predicciones %>% 
  bind_cols(data_test_prep %>% dplyr::select(yield)) %>%
  dplyr::mutate(modelo = "lm",
                yield = exp(yield),
                .pred = exp(.pred))

summary(predicciones_rf)  
     .pred          yield         modelo         
 Min.   :2452   Min.   :2380   Length:236        
 1st Qu.:5084   1st Qu.:5118   Class :character  
 Median :6068   Median :6104   Mode  :character  
 Mean   :6061   Mean   :6060                     
 3rd Qu.:7045   3rd Qu.:7035                     
 Max.   :8979   Max.   :8969                     

Testeo interno

Resumen de las métricas de predicción

error_test_lm  <- multi_met(
  data     = predicciones_lm,
  truth    = yield,
  estimate = .pred,
  na_rm    = TRUE
) %>%
  dplyr::mutate(
    modelo = "lm"
  )
error_test_lm
# A tibble: 4 × 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard     140.    lm    
2 rsq     standard       0.990 lm    
3 mape    standard       1.89  lm    
4 mae     standard     108.    lm    

Testeo interno

Discusión final

  • Éxito = Conocimiento del Negocio + Conocimiento estadístico + Capacidad de resolver problemas.
  • Cada modelo necesita de una receta de preprocesamiento adecuada.
  • Las predicciones deben estar en el mismo dominio que la variable observada.
  • Muchas variables => maldición de la dimensionalidad => sobreajuste.
  • No existe una sola forma de resolver un problema.

¿Preguntas?